Prediction of equilibrium Li isotope fractionation between minerals and aqueous 
solutions at high P and T\ an efficient ab initio approach 

Piotr M. Kowalski and Sandro Jahn 

GFZ German Research Centre for Geosciences, Telegrafenberg, 14473 Potsdam, Germany 



o 

(N 



(N 

o 

o 

Ph. 



(N 

en 

in 

o 



X 



Abstract 

The mass-dependent equilibrium stable isotope fractionation between different materials is an important geochemical 
process. Here we present an efficient method to compute the isotope fractionation between complex minerals and 
fluids at high pressure, P, and temperature, T, representative for the Earth's crust and mantle. The method is tested by 
computation of the equilibrium fractionation of lithium isotopes between aqueous fluids and various Li bearing miner- 
als such as staurolite, spodumene and mica. We are able to correctly predict the direction of the isotope fractionation 
as observed in the experiments. On the quantitative level the computed fractionation factors agree within 1.0 %c with 
the experimental values indicating predictive power of ab initio methods. We show that with ab initio methods we are 
able to investigate the underlying mechanisms driving the equilibrium isotope fractionation process, such as coordi- 
nation of the fractionating elements, their bond strengths to the neighboring atoms, compression of fluids and thermal 
expansion of solids. This gives valuable insight into the processes governing the isotope fractionation mechanisms on 
the atomic scale. The method is applicable to any state and does not require different treatment of crystals and fluids. 
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1. Introduction 

The fractionation of stable isotopes between var- 
ious materials is of importance in geoscience, as 
the variation in isotope content provides valuable in- 
formation on processes and interaction between at- 
mosphere, biosphere, geosphere and hydrosphere. 
Although there is substantial analytical work per- 
formed in this area, reliable computational meth- 
ods to predict isotope fractionation factors have been 
available only recently, proving that they can con- 
tribute towards understanding geochemical mecha- 
nisms responsible for production of isotope signatures 
JDriesner, 1997; Yamaiietal. 2001; Schauble, 2004 
i Domagal-Goldman et all l2008t Iffijl & Schauble, 2008 
Meheut et al.','2009VSchauble et al.','2009; Zeebe, 2009 
Hilletal., 2010; Rustad etal., 2010a; Rustad et al, 
l20 1 Obt IZeebi l2009[ l20 1 Ol) . 

Ab initio calculations of equilibrium isotope frac- 
tionation between minerals have received considerable 
attention recently. Previous studies, however, were 
mostly limited to simple crystals containing just a few 
atoms in the unit cell such as quar t z, kaolinite or carbon - 
ate minerals (iMeheut et all 120071 iRustad et all l2010ah . 



as the methods used require considerable computational 
resources. Only very recently, the calculations have 
been extended to more complex crystalline s olids con- 
taining up to 80 atoms in the unit cell by [Schauble 
(120111) . There are different approaches used in the com- 
putation of the mass-dependent stable isotope equilib- 
rium fractionation factors of minerals, but all methods 
require knowledge of the vibrational spectrum of the 
considered sys tem, which is usually computed using ab 
initio methods. iMeheut et all (120071) performed full nor- 
mal mode analysis of the solid phases accounting for 
the phonon dispersion in reciprocal space. Because of 
the huge computational requirements, this method, al- 
though correct, can be only applied to the computation 
of stable isotope fractionation between simple phases. 
On the other hand, in order to derive the frequencies 
requ ired for the computa tion of the fractionation fac- 



tors 



Rustad et al.l (120 10 a) approximated the solids by 



small clusters and treated them as large molecules. This 
approach is based on well estabhshed theories of sta- 
ble isotope f ractionation ( Bigeleisen & Mayeij Il947 : 



Kieffer et al.L 119821: IChackoet a l.. 2001) showing fliat 



the major contribution to the mass-dependent fraction- 
ation comes from the local vibrational motion of the 
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Table 1: Lattice parameters of the investigated Li-bearing silicates. Naioms is the number of atoms in the modeled supercell. The units are A and 
degrees. 

staurolite spodumene mica IM mica 2M1 mica 2M2 mica 3T 



a 
b 
c 
a 
P 

y 

ref. 

N„tn, 



7.848 

16.580 

5.641 

90 

90 

90 

1 

81 



9.463 
8.392 
10.436 

90 
110.15 

90 

2 

80 



5.20 
9.01 
10.09 

90 
99.38 

90 

3 

44 



5.209 
9.053 
20.053 

90 
95.74 

90 



9.04 

5.22 

20.2100 

90 
99.58 

90 



5.200 

5.200 

29.760 

90 

90 

120 

5 

66 



References: IComodi et alJ (l2002l) . fcameron et al.1 (ll973h . ^Sartoril (ll976h . 1Sartori et all (ll973h . terownl (Il978h 



fractio nating element. In line with this finding lSchauble 
(1201 ih has found that considering the phonon spectrum 
on a single phonon wave vector only is sufficient for 
modeling of ^^Mg/^'^Mg isotope fractionation between 
magnesium bearing crystal phases. 

As fluid-rock interactions are a major cause that 
alter the isotopic signature of a mineral in a rock, 
understanding the equilibrium isotope fractionation 
processes between minerals and aqueous fluids is of 
great importance in petrology. Although there has 
been considerable work on stable isotope fractiona- 
tion between various minerals and the computational 
techniques are well established, the question of 
treating fluids, namely aqueous solutions remains 
open. Most of the ab initio calculations of iso- 
tope fractionation in flui ds use t he cluster approach 
JDomagal-Goldman et jL[__|2008'; 'Hill & SchaubT 



2008: Zeebe. 2009; Rustadetal, 2010b; Rustadetal. 



201 Oat lHiUetal.1 l2010t lYamaii et all I2OOII: IZeeb" 
2010V in which the considered species (ions or molec- 
ular complexes such as Fe, Mg, H3BO3) are surrounded 
by a hydration shell and the whole structure is relaxed 
assuming T = K. This approach is based on the 
computation of static atomic configurations and is valid 
at low temperatures only. In case of Li in aqueous 
solution at high temperatures {T ~ 1000 K), frequent 
exchange between particles of the hydration shell 
surrounding Li cation with the fluid is observed on time 
scales as short as picoseconds (10"'^ 
(l2009l) V Distribution of cation coordination 



Jahn & Wunder 



and 



cation-O bond lengths, effects that are expected to 
aff'ect the isotope fractionation (Bigeleisen & Maver, 



1947) 



2009 



also change with pressure (Jahn & Wunder, 
Wunder etal., 2011). These features are difficult 



to account for by using the cluster approximation for a 
compressible fluid at high temperature. The impact of 
the dynamical behavior of particles and compressibility 



of fluid must be investigated in order to properly 
compute the isotope fractionation in aqueous fluids. 
The only recent ab initio work that accounts for the 
dynamical efi'ects on the isotope fractionation in fluid 
is by Rustad & Bvlaska (2007) who considered boron 
isotope fractionation between B(OH)3 and B(OH)4 in 
aqueous solution. They performed ab initio molecular 
dynamics simulations of this system and tried to use the 
vibrational density of states derived through the Fourier 
transform of the velocity auto-correlation function 
as an input for the calculation of the "B/'°B isotope 
fractionation coefficient. The resulting fractionation 
factor a = 0.86 is much lower than the experimental 
value a = 1.028. Interestingly, the discrepancy between 
experiment and theory is cured by quenching the 
selected configurations along the molecular dynamics 
trajectory and computing the harmonic frequencies. 
The fractionation factor derived using these frequencies 
exactly reproduces the experimental value. 

In this contribution we present an efficient approach 
to the computational prediction of equilibrium isotope 
fractionation between complex minerals and fluids at 
high P and T. Both solids and fluids are treated as 
extended systems by application of periodic bound- 
ary conditions in all three spacial directions. We will 
demonstrate that at T > 600 K the fractionation factor 
can be computed by considering the force constants act- 
ing on the fractionating element only. Both solid and 
fluid supercells should be big enough to avoid signif- 
icant interaction between atoms and their periodic im- 
ages. In our investigation we use cells at least 5 A wide 
in each spacial dimension. A representative distribution 
of relevant coordination environments in the fluid struc- 
ture is obtained by perfor ming Car-Parrinello mo lecu- 
lar dynamics simulations (^Car & Parr inelloill985b . For 
the calculation of the fluid fractionation factors, several 
random snapshots from this simulation are chosen. The 
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Figure 1: 
spectram 



The/? factors for staurolite and spodumene. The hnes represent the results using Eq. [5](solid lines) and Eq. [T]with the full frequency 
(dotted lines). 



force constants acting on the fractionating element and 
the resulting fractionation factors are then obtained for 
each configuration and the fractionation factor for the 
considered element in the fluid is computed as an aver- 
age over the whole set of geometries. 

As a test case for our approach, we have computed the 
fractionation factors between Li bearing aqueous fluids 
and three minerals, mica, staurolite and spodumene. For 
these system s, recent exper i mental data a re available for 
comparison ( IWunder et all l2006l 120071) . Furthermore, 
lithium as one of the lightest elements with two stable 
isotopes produces strong isotope signatures. It strongly 
fractionates into aqueous fluids during fluid-rock inter- 
action processes and is use d as a tracer of mass trans- 
fer in the subduction cycle JWunder et al.l 120061) . The 
two stable isotopes, ^Li and ^Li, have respective abun- 
dances of 92.5% and 7.5%. The large mass difference 
of 7.016003/6.015121 = 16.6% results in a prominent 
fractionation of at least a few %o even at high temper- 
atures T ~ 1000 K. The experimental data on Li iso- 
topes indicate a significant influence of the Li coordi- 
nation and the Li-O bond length on the fractionation 
of Li isotopes. The heavier isotope preferentially occu- 
pies the lower coordinated sites and phases with shorter 



bond distance (IWunder et all 120111) . w hich is expected 



also f rom the theoretical point of view JSchauble et al. 
20091) . We will show that the application of ab ini- 



tio methods to Li-bearing silicates and fluids provides 
unique insight into the mechanisms driving equilibrium 
Li-isotope fractionation on the atomic scale. 



2. Theoretical model 

The mass-dependent equilibrium isotope fractiona- 
tion is driven by the change in the molecular and crys- 
talline vibration frequencies resulting from the different 
mass of the isotopes. The fractionation between species 
and an ideal atomic gas is called the /3 factor or the re- 
duced partition function ratio (RPFR) and in the har- 
monic approximation is given by the formula: 



nM,. 



exp 



(m; - u*) 1 - exp(-M,) 
2 1 - exp(-M*)' 



(1) 



where u = haji/ksT, h = h/ln is the reduced Planck 
constant, w,- the vibrational frequency of the i-th de- 
gree of freedom, ks is the Boltzmann constant, Ndof 
is the number of degrees of freedom, which for the A^ 
being the number of atoms in the considered system 
(molecule, mineral or fluid) is equal to 3A^ - 5 for a 
diatomic molecule, 3A^ - 6 for multiatomic molecules 
and 3A^ for crystals, and a star symbol marks the heav- 
ier isotope. Despite requiring only the knowledge of 
the vibrational frequency spectrum, the above formula 
accounts also for th e translational and r otational mo- 
tions of a molecule (IChacko et al.1 l200lh . Because of 
the Redlich-Teller product rule, equation[T]is also valid 
for minerals (but with the product runnin g to 3A^), if the 
crysta l is represented as a big molecule dChacko et al. , 
2001i) . The fractionation factor between two substances 
A and B, qa-b is computed as the ratio of the relevant/? 
factors, which for (J5 - I) ~ 10"^ is well approximated 



Table 2: The fi factors for mica (columns 2-5) and fractionation factors between mica and spodumene (last column) computed at T = 650 K for 
various mica polvtvpes and Li substitution sites. All values are given in %o. The measured value is that of Wunder et alj i2007t). 


Mineral 


Lil 


Li2 


Li3 


Average 


A Limc.-spd. 


IM 


13.9 


14.9 


- 


14.6 


+4.7±0.9 


occupation 

2M1 


0.3 
13.9 


0.97 
13.6 


- 


13.7 


+3.8±0.9 


occupation 

2M2 


0.38 
13.8 


0.92 

13.4 


- 


13.6 


+3.7±0.9 


occupation 
3T 


0.37 
12.2 


0.95 
14.9 


12.1 


13.6 


+3.7±0.9 


occupation 
exp. 


0.7 


0.89 


0.14 




+2.5±1.0 


by its differences: 






tions can be 


significant as the full normal mode analysis 



ClA-B = PaIPb - ^A-B = Pa -Pb- 



(2) 



The calculation of the y6 factor requires only the knowl- 
edge of the vibrational properties of the considered sys- 
tem computed for the two different isotopes. How- 
ever, computation of the whole vibrational spectra of 
complex, multiparticle minerals or fluids requires sub- 
stantial computational resources and is currently lim- 
ited to systems containing a few dozens of atoms. 
Any approach that would allow for a substantial re- 
duction of computational time and computationally ef- 
ficient treatment of complex systems is highly desired. 
iBJg eleisen & Maver (1947) have shown that in case of 
M < 2 the isotope fractionation can be computed from 
the knowledge of the force constants acting on the atom 
of interest. The yS factor (Eq. [TJ can then be approxi- 
mated by: 



p^i+Y. 



^'^«? 



24 



= 1-H 



Am 



mm* lAklT^ . , 



Yj^iO) 



where A,- are the force constants acting on the isotopic 
atom in the three perpendicular spacial directions (x, y 
and z). Am - m* - m and m is the mass of the frac- 
tionating element. For clarity we will call the formula[3] 
the single atom approximation through the paper The 
validity criterion, u - hoj/ksT < 2, restricts the usage 
of the formula to frequencies (w [cm ' ] < 1.39 r[K]. As 
it is rare that ai » 1000 cm"' (with the exception of 
the vibrations involving hydrogen), the formula is usu- 
ally valid for high temperatures T > 600 K. In the case 
of Li, o) < 600 cm ' and the formula is valid down to 
T ~ 450 K. This gives us the opportunity to simplify 
the calculations by considering the force constants act- 
ing on one atom of interest instead of all atoms consti- 
tuting the considered system. For large systems con- 
taining hundreds of atoms the speed up in the calcula- 



of an A^-atoms system requires A^ times more computa- 
tions than computing a single atom. For instance for a 
system containing 100 atoms the calculations using the 
single atom approximation are 100 times faster We will 
show that the computation of isotope fractionation fac- 
tors from the knowledge of the force constants acting 
upon the element of interest allows for efficient calcu- 
lation of Li isotope fractionation between complex sil- 
icates, such as spodumene, Li-micas and Li-staurolite, 
and aqueous solutions. 

One important aspect of the method is its usage for 
the calculation of isotope fractionation in crystals. In 
principle in order to compute the yS factors for crystals 
one has to account for dispersion. In a solid the phonon 
frequencies are identified by a (7-vector in a reciprocal- 
space, which requires extension of the product in Eq. 
[1] beyond the number of atoms and addin g multiplica 



tion ov er the ^-vector grid (see Eq. 16 of iMeheut et al. 
(l2007b ). However, considering the ^^Mg/^'^Mg frac 



tionation in Mg-bearing minerals ISchaubld (120111) has 
shown recently that for minerals of multiatomic struc- 
ture (A^ > 20) considering a single phonon wave-vector 
is sufficient for getting very accurate yS factors even at 
r = 300K (eiTor of 0.1 %c). At T = 1000 K the error 
is negligible and in the order of 0.01 %c. This finding 
and the computation of yS factors considering the sin- 
gle atom approximation reduce the computational load 
required to compute the fractionation factors to calcu- 
lation of only the force constants acting upon fraction- 
ating element. This allows for computer-aided investi- 
gation of isotope fractionation in complex minerals and 
fluids containing hundreds of atoms. 

3. Computational approach 

The calculations of p factors of crystals and aque- 
ous solutions were performed by applying density func- 
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Figure 2: The/J factor for various polytypes of mica. Tlie lines represent the results for isotope substituted on different Li sites: Lil (solid lines), 
Li2 (dotted lines) and Li3 (dashed line). 



tional theory (DFT) methods, which are currently the 
most efficient methods allowing for treating extended 
many particle systems quantum-mecha nically. For that 
purpo se we used the CPMD code (iMarx & Hutter . 
20001) . which is especially suited for ab initio simula- 
tions of fluids. In order to reach consis tency with previ- 
ous w ork on Li-bearing aqueous fluids jjahn & Wunder , 
20091) . we used the BLYP exchange-c orrelation func- 



tional jBeckel Il988c iLee et al.L Il988h . a plane wave 



basis set and an energy cut-ofF of 70Ryd for geom- 
etry relaxations and molecular dynamics simulations 
and of 140 Ryd for computation of vibrational frequen- 
cies and force constants. The much higher cut-off 
used for derivation of the vibrational frequencies and 
force constants was essential to obtain the converged /3 
factors. Norm-conserving Goedecker pseudopotentials 
were apphed for the d escription of the core electrons 



JGoedecker et al.Lll996) . For both crystalline solids and 



aqueous solutions, periodic boundary conditions were 
applied. The solids were represented by large cells con- 
taining at least 40 atoms. The number of atoms used 



in the crystal calculations together with the lattice pa- 
rameters of modeled crystals are summarized in table[T] 
The lattice constants used in our calculations resemble 
those determined bv I Wunder eT al. (2006, 200"^. The 
atomic positions of the crystal structure were relaxed to 
the equilibrium positions to minimize the forces acting 
on the atoms. The aqueous solution was represented 
by a periodically repeated box containing up to 64 wa- 
ter molecules and one Li atom. The Li^ cation in the 
fluid was charge balanced by an F" anion. The pres- 
sure and temperature conditions w ere chosen to be close 
to the experimental conditions of lWunderet alj (120061 
20071) . The pressure of aqueous solution for a given 



temperature and volum e was calculated according to 
the equation of state of IWa^ner & PrussI (12002) . The 
ab initio molecular dynamics simulations (AIMD) were 
preformed for fixed temperature and yolurn e using Car- 
Parrinello scheme dCar & Parrinellol Il985h . The tem- 
perature during each ru n was controlled by a Nose- 
Hoov er chain thermostat (INose & Kleiru.ll983tlHoover . 
1985i) . For each T -V conditions at least 10 ps long tra- 
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Figure 3: The structures of [Li(H20)n]^ clusters. 



jectories have been generated with an integration step of 
0.12 fs. The sum of the force constants needed for com- 
putation of thCyS factors from equation |3] was computed 
using finite displacement scheme by fixing the positions 
of all the atoms except the fractionating element. The 
full normal mode analyses were performed using the 
same method, but displacing all the atoms constituting 
the considered system. The frequencies were obtained 
through the diagonalization of the full dynamical ma- 
trix (ISch auble. 2004) as implemented in CPMD code. 
In case of solids the atomic structures taken for compu- 
tations of j6 factors were those obtained after relaxation 
of atomic positions to minimize the forces for given lat- 
tice constants. For fluids the/? factors were computed on 
the ionic configuration snapshots extracted uniformly in 
0.1 ps intervals along the 10 ps long molecular dynam- 
ics trajectories. The calculations were performed with 
the positions of water molecules fixed to the molecular 
dynamics configurations and the Li cation was relaxed 
to the equilibrium position. The effect of the continuous 
medium on the derived fractionation factors was studied 
by additional computations of Li(H20)|^ isolated clus- 
ters. For that purpose we used a large, isolated simula- 
tion box of the length of 16 A, forcing the charge density 
to be zero at the boundary, as implemented in CPMD 
code. 

The error in the computed value of the /3 - I and A 
fractionation factors we estimate from an average error 
of vibrati onal frequencies computed using chosen DFT 
method. Fi nlev & StephensI (Il995b : iMenconi & Tozer 
(12002) estimated the errors made in calculations of vi- 
brational frequencies of small molecules using differ- 
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Figure 4: The /3 factors for [Li(H20)„]+ clusters. The lines represent 
the results for n=3,4,5 and 6 (from top to bottom) of this work (solid 
lines) and using frequencies computed in Yamaji et al. 1 2001) (dashed 
lines). The results for n=3 and n=4 are nearly identical and hardly 
resolved in the figure. 



ent DFT functionals. According to these works BLYP 
functional systematically overestimates the frequencies 
by ~ 3.5 % with the deviation from the mean value of 
~ 1%. Therefore, we expect that using BLYP func- 
tional the j3 - 1 and A values are systematically over- 
estimated by 7 % and that in addition there is 2 % error 
in derived yS - 1 factors. We notice that in order to cor- 
rec t for the system atic errors some authors (for exam- 
ple lSchaubldOOl II) ) scale the DFT vibrational frequen- 
cies usually by a frequency independent scaling factor, 
which could be derived from the match to the experi- 
mental measurements. We decided not to use such a 
scaling as we intent to test the ability of DFT methods 
to predict the stable isotope fractionation factors from 
first principles without introducing free parameters, or 
making constraints to the experimental data. 



4. Results and discussion 



4.7. Solids 



4.1.1. Representation of the silicates 

The lattice parameters of the modeled crystalline 
solids are the experimental values found in the liter- 
ature. For staurolit e, the refined crystal structure of 
Comodi et al. ( 2002h w as used. As in the experiment 
of IWunder et al.l (120071) Mg- staurolite was used instead 
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Figure 5: The dependence of the /3 factors of [Li(H20)„]+ clusters on 
th e Li-0 d istance in the clusters derived using vibrational frequencies 
of lYamaii et al. (2001) (circles) and ours (squares). The coordination 
runs from 3 (left) to 6 (right). The dotted lines connecting the data 
points are added to visualize the trend. 



of Fe-staurolite, in order to reproduce closely the ex- 
perimental conditions we replaced all the Fe atoms in 
the modeled structure with Mg atoms. In staurolite 
Li is a trace specie s. Following the assignment of 
Wunder et alJ (l2007h we assumed that it occupies one 
of the T2 sites, i.e. the 4-fold coordinated site occu- 
pied by Mg atoms, and that there is only one substitution 
site. The constructed model contains a single unit cell of 
chemical composition Al'j';^(LiiMg3)t'^lSi^'*'045(OH)3, 
where in square brackets we denote the coordination 
number. The chosen composition and lattice parameters 
cl osely resemble s the on es determined for Mg- staurolite 
bv lWunder et"ai] (l2007b . 

Spodumene is the simplest cr ystal investigated here . 
The modeled structure is that of ICameron et al.l (119731) . 
The chemical composition of the unit cell used in the 



investigation is (Li8Alg)'^''lSijg048, which is exactly the 
chemical compositions of spodumene sy nthesized and 
used in the isotopic measurements by IWunder et al. 



(I2006h . 

Comparing 
Li-bearing 



with staurolite and spodumene the 
mica obtained in the experiments by 



Wunder et al.l (120071) is a complex silicate system con- 
taining different polytypes with relative abundances 
varying significantly between different samples (see ta- 
ble 3 of Wunder et al. (2007)) . Foll owing structure 



gation we consider four mica polytypes: IM, 2M1, 
2M2 and 3T. The structural parameters and litera- 
ture sources are given in table [1] In order to model 



bleJTL 
IWund 



determination of Wunder et al. (2007) in our investi 



the minerals synthesized in IWunder et al.l (120071) ex- 
periment we represent the different mica polytypes 
by the supercells of the following chemical compo- 
sitions: K2(Li4Al2)[*lSif'*^02o(OH)4 for IM mica, by 
K4(Li8Al4)[*]Sif'J,l04()(0H)8 for 2M1 and 2M2 micas, 
and K3(Li6Al3)[''lSi'i2'03()(OH)6 for 3T mica. 

4.1.2. Li isotope fractionation between minerals 

StauroUte and spodumene have a single Li occupation 
site. In staurolite, Li substitutes for Mg in a four-fold 
coordinated site, while in spodumene Li occupies the 
six-fold coordinated M2 site of pyroxenes. In both cases 
Li is bounded to oxygen atoms only. The computed fi 
factors for both silicates are given in figure [T] We give 
the results of two sets of calculations: (1) considering 
force constants acting upon Li atom only using Eq. |3] 
and (2) performing full normal mode analysis, i.e. com- 
puting full phonon spectrum at the gamma point and us- 
ing Eq. [T] This provides an explicit test of the single 
atom approximation outlined in section |2] The yS factors 
derived using both methods are essentially identical and 
only deviate slightly at low temperat ures, which is ex- 
pected. Wunder et al. (2007) and Wunder et all (l2006h 
measured the Li isotope fractionation between these two 
minerals and the aqueous solution. According to their 
measurements the fractionation between staurolite and 
spodumene is 2.7 + LO %c at 1200K and 3.7 + LO %o at 
1000 K. The calculated values, which are given by the 
differences between /? factors at the considered temper- 
atures are A^Lisn-.-spd. = Atr. -Apd. - 3.7 + 0.5 %o and 
4.6 + 0.5 %c respectively and therefore in good agree- 
ment with the experiment. We will show that because 
of thermal expansion effect and different experimental 
pressures (3.5 GPa with staurolite and 2.0 GPa in exper- 
iments with spodumene), the computed A^Lijtr.-spd. is 
overestimated by 1 . 1 %o, bringing the prediction to even 
better agreement with the measured data. 

The case of mica is more complex as it contains dif- 
ferent p olytypes and L i substitution sites. In the exper- 
iment of IWunder et al.l (2007) the measured mica sam- 
ples contained various combinations of IM, 2M1, 2M2 
and 3T polytypes. To account for that we have com- 
puted the p factors for all the outlined polytypes and Li 
substitutions sites. The results are given in figure|2] It is 
clearly visible that both the polytype and Li substitution 
sites impact slightly the value of computed y6 factors. 
This is because the different structural environments re- 
sult in slightly different Li-0 bond lengths, although the 



Table 3: The distribution of coordination number of Li in aqueous 
solutions computed by Jahn & Wunder (2009). 
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0.28 
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0.52 
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0.00 


0.75 


0.31 


0.53 


0.15 


0.01 


1.2 
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0.54 


0.21 
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0.24 


0.01 
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Li coordination is the same in all cases. The largest dif- 
ference is visible in case of 3T polytype, where j3 factor 
computed for Li2 site is higher than for the other Li sites 
and polytypes. This is because even after atomic relax- 
ation this particular site exhibits the shortest Li-0 bonds 
with the strongest Li-O bond shorter by ~ 0.05 - 0.1 A 



compa ring to other Li sites and polytypes. IWunder et al. 
(l2007h showed that at approximately T - 650 K the 



fractionation between mica and spodumene minerals is 
2.5 ± 1 %c. The results of our calculations for that tem- 
perature are reported in table |2] Here we derived the fi 
factors for the different mica polytypes by taking the sta- 
tistical average over the /? factors computed for each Li 
site. The contribution of each site is weighted according 
to the occupation of the particular site by Li atom, which 
is also given in table|2] The calculations predict the cor- 
rect fractionation direction, i.e. A^Limca-spd. > and 
the experimental fractionation factor within uncertain- 
ties of the calculations (which we estimate at ~ 0.9 %o 
at considered temperature) but slightly overestimate the 
measured value. We will show later in the discussion of 
the fractionation between solids and fluid that account- 
ing for thermal expansion of the crystals the reported 
computed values decrease by ~ 0.3 %c further improv- 
ing the agreement with the experiment. 

4.2. Fluid 

4.2.1. Cluster approach 

In most of the recent work on ab initio computa- 
tion of the stable isotope fractionation in aqueous so- 
lutions the isolated cluster approach is used in which 
a considered species is surrounded by the hydration 
shell and the whole structure is optimized assuming 
T = K dYamaii et all bOOll iDomagal-Goldman et al. 1 
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Figure 6: The pressure dependence of the fi factor for Li in the 
fluid computed based on the cluster approach using the vibrational 
frequencies of Yarnaii et aL (2001) (circles), the full frequency spec- 
trum computed for c luster s in this work (diamonds) and using the 
iBigeleisen & Mavea 1 19471) approximation (their Eq. (21)) tog ether 
with Li-0 symmetric stretching frequencies of Yamaji et al. 1 2 0011) 
(squares). The bars represent the yS factors computed along the ab ini- 
tio molecular dynamics trajectories and their width represent the un- 
certainties in computed values. The dotted lines connecting the data 
points are added to visualize the trend. 



high temperatures and pressures the hydration shell sur- 
rounding lithium i on is not static but exhi bits strong dy- 
namical character (IJahn & \yunden.l2009h and co mpres- 



sion impacts its structure d Wunder et all 120111) . The 



important questions are on the impact of these effects 
on the equilibrium isotope fractionation and how well 
these effects can be described with the widely used clus- 
ter approach. In order to address these problems we per- 
formed set of calculations involving [Li(H20)n]^ clus- 
ters. The clusters used in the investig ation are illustrated 
in figure |3] Following the work of Yamaii et al.l (1200 lb 
we computed the p factors for isolated [Li(H20)n]"^ 
clusters for n = 3,4,5,6, relaxing the structures to 
equilibrium positions and computing the full frequency 
spectra. The spectra were then used to compute /? 
factors according to Eq. [1] In the same way we 
als o computed the 13 fac tors using frequencies derived 
by lYamaii et al.l (1200 ih obtained with the restricted 
Hartree-Fock method (RHF). Both results are given in 
fi gure m The fi factor s computed with the frequencies 
of lYamaii et al.l (1200 lb are higher than the ones derived 
with DFT frequencies except in the n - 5 case, for 
which both calculations predict the same values. This 
may be related to different cluster structures used in the 
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Figure 7: The dependence of the fi factor of fluid on the size of the 
simulation cell. The thick and thin dashed hues represent the aver- 
age value and the uncertainty limits of /J factor computed on system 
containing (slUiO molecules. 



calculations as pos itions of hydrogen at oms are not pro- 
vided in details bv lYamaiiet alj (1200 lb . An interesting 



observation is illustrated in figure |5] where the fi fac- 
tor is plotted as a function of Li-O bond length. With 
increasing n the Li-0 bond length increases, as the wa- 
ter shell containing more water molecules has to relax 
outwards creating more space for additional molecules. 
The increase in the bond length results in a decrease of 
the jS factor. This has an important implication on the 
pressure dependence of the fi factors derived using the 
cluster approach. 

Having both results for clusters we attempted to in- 
vestigate the pressure effects on the fi factors. We do 
that by averaging the fi factors over the statistical distri- 
bution of [Li(H20)n]^ comple xes in aqueous solution , 
which is pressure dependent. Ijahn & Wunden (120091) 
have shown that in the pressure range from 1 to 6 GPa, 
the Li coordination by oxygens increases smoothly from 
preferentially four-fold to five-fold coordination. At 
2 GPa, which corresponds to the experimental condi- 
tions of Wunder et al. (2006, 2007) , the mean Li co- 
ordination i s about 4.2. We took the probability dis- 
tribution of Ijahn & Wundea (l2009h . which is given in 
table [3] and derived the pressure-dependent fi factors 
as a statistically weighted average of the fi factors de- 
rived for [Li(H20)n]^ clusters. The results are given 
in figure |6] Both results derived on the two fi fac- 
tor estimations predict decrease of the Li isotope frac- 
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Figure 8: The /? factors for Li^ aqueous solution at 1 .9 GPa obtained 
by [Li(H20)n]^ cluster approach and AIMD simulations. The lines 
represent the results of this work (solid l ine) and the fi factors ob- 
tained using frequencies of I Yamaji et al J ^2001 ) (dotted line). The 
dashed line represents the fit to the /3 factors given by bars, which 
indicate the uncertainty in the computed values, and computed at dif- 
ferent temperatures as an average over the AIMD trajectories. 



tionation with increase in pressure. This is because at 
higher pressure the more coordinated and with longer 
Li-O bond lengths [Li(H20)„]^ structures are preferred, 
which results in lower/? factors. This finding is counter 
intuitive, as one should expect that the compression 
of the fluid should lead to the shortening of the Li-O 
bonds, elevated vibrations and resulting higher fi fac- 
tors. In figure |6] we also give the estimation of fi fac- 
tors computed from the knowledge of the Li-O totally 
symmetric stretching frequencies using rough approx- 

3^ (their Eg. 21} 



imation of iBigeleisen & Maverl d 1 94 



(1200 Ih 



with the relevant frequencies of Yamaii et al 

an d the [Li(H20) i 1^ clu sters pr obability distribution 



of Ijahn & Wundej (l2009t) . In the IBigeleisen & Maver 
(1194 7) approximation the fi factor is proportional to the 
square of the totally symmetric stretching frequency, 
Vs, and the cluster size, i.e. fi ~ v^n. As with in- 
creasing the cluster size, Vj decreases by ~ 20%, the 
largest effect on the isotope fractio nation computed us- 
ing the lBigeleisen & Maverl (119471) method comes from 
the coordination (cluster size). The resulting pressure- 
dependent fi factor shows the desired tendency. It in- 
creases with the size of the cluster, which causes the 
increase in pressure as is seen in figure|6l We will show 
that the simulation of continuous media is required for 
proper investigation of the effect of the compression and 
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Figure 9: The average Li-O distance between Li and the three closest 
O atoms in aqueous fluid as a function of pressure. 



to obtain realistic isotope fractionation signature of high 
P fluids. 

4.2.2. Molecular dynamics approach 

In order to fully account for the pressure effects, spa- 
cial continuity of the fluid and its dynamical motion 
we produced lOps long molecular dynamics trajecto- 
ries of systems consisting of 64 H2O molecules and 
one Li ion for different T = lOOOK, 800 K and 600 K 
and pressure of 1 .9 GPa, whic h closely resembles the 
experimental conditions of Wunder et alj (2006, 20071) . 
The corresponding simulation box length is 12.17 A at 
T - 1000 K. We note that the thermal eff'ects on the 
pressure will require to use a supercell of ~ 1 % smaller 
box length for T - 600 K, a small eff'ect which we omit- 
ted. As the oxida tion state of Li in the aqueous solution 
is -I- 1 , following Jahn & Wunden (12009) we added a F 
atom to the system as a charge compensator An inter- 
esting question is on the impact of the system size on the 
computed j6 factors. In order to investigate this problem 
we computed lOps length trajectories also for simula- 
tion boxes containing 8, 16 and 32 H2O molecules for 
T - 1000 K and pressure of 1.9 GPa. The resulting /? 
factors are given in figure |7] Within the accuracy of 
the calculation the fi factor is system size independent 
and in principle small systems containing 8 H2O atoms 
can be used in the investigation. This substantially re- 
duces the required computational time. As the current 
implementations of plane-wave DFT methods scale as 
A^^ -A^^, with A^ being the number of particles in the sys- 



tem (number of electrons), the computation time gained 
reducing the number of particles in the computational 
box could be significant. In our calculations switching 
from a system containing 64 water molecules to 8 the 
gain is a factor of ~85. Nevertheless for our calcula- 
tions we used the simulation box containing 64 water 
molecules. In order to obtain the temperature dependent 
fi factor at P = 1.9 GPa we fitted by the least squares 
procedure the formula j6 - l+AjT^ to the /J factors cal- 
culated at the three temperatures. The parameter of the 
fit is A = 6. 112 ■ lO""* for temperature expressed in units 
of 10^ K. The resulting fi factor as a function of temper- 
ature at P = 1.9 GPa is given in figure [8] together with 
the already discussed predictions using the cluster ap- 
proach. Interestingly, the molecular dynamics fi factor 
is in good agreement with the value obtai ned by using 
clusters approach with Yamaii et al.l (12001 ) frequencies. 
The difference between our cluster and MD calculations 
is also moderate, 0.6 %o at 1000 K and 1 %o at 800 K. 
However, the agreement between both types of calcu- 
lations is only reached at lower pressures (P < 2 GPa), 
which will be discussed in the next paragraph. In or- 
der to check the validity of the single atom approxima- 
tion outlined in section |2] for fluids we computed/? fac- 
tors with the full frequency spectra obtained for selected 
configurations. We found only negligible deviation of 
the resulted j6 factors from the ones derived considering 
force constants acting on the Li atom only. 

In the [Li(H20)„]^ cluster calculations we obtained 
an unexpected result indicating that/? factor should de- 
crease with pressure, which we found counter intuitive. 
In order to investigate the pressure impact on the fi fac- 
tor accounting for the continuity on the medium and 
its pressure-driven compression we computed the frac- 
tionation factors at r = lOOOK and different pres- 
sures on a system containing 8 H2O molecules. In each 
case 10 ps long trajectories were generated and fi fac- 
tors were computed on a set of atomic configurations 
extracted uniformly along the trajectories. The result 
is given in figure |6l We clearly see that for pressures 
P > 2 GPa, as the effect of compression, the fi factor 
increases monotonically with increasing pressure. The 
reason for that is the small decrease in the mean Li- 
O bond length (measured as the average distance be- 
tween Li and the three closest neighbors) with increas- 
ing P, which is opposite to the result using clusters ap- 
proach, and the coordination, as shown in figure|9] This 
finding is in line with results of Wunder et al. (201 lb , 
who found that the mean Li-O distance increases for the 
pressures up to 1 GPa and remains constant at higher 
pressures. This explains why the computations using 
cluster approach, in which the increase of the Li-O bond 
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Figure 10: The solid-fluid isotope fractionation between Li-bearing minerals and aqueous solution. The lines represent the computed values and 
the shadowed areas reflect the uncertainties computed assuming the computational error evaluation described in section[3] IM, 2M1, 2M2 and 3T 
lines indicate results obtained for different mica polytypes. The points are the measured values of Wunder et al. (2006, 2007). Left and right panels 
represent the result without and with thermal expansion correction discussed in the text. The thermal expansion coirection is made by a constant 
negative shifts of the solid-fluid fractionation factors of -0.6 %o for staurolite, -0.7 %o for mica and -0.4 %o for spodumene. The pressure assumed 
for the computations is 2 GPa . The computed values for staurolite given in right panel are shifted by additional -0.5 %o, the correction due to the 
higher experimental pressure P = 3.5 GPa, as discussed in the text. 



lengths with the increase in the cluster size, and there- 
fore pressure, is also observed, produce good pressure 
dependence of yS factor at low pressures, as illustrated in 
figure |6] On the other hand this clearly shows that an 
isolated cluster is not a good representation of high P 
and T fluid and can not be used for the computation of 
theyS factors in fluids at extreme conditions. Continuity 
and compressibility of the fluid have to be considered in 
order to obtain realistic results. 

Although most of the experimental results to which 
we refer in this paper were performed at lower pres- 
sures (2-3.5 GPa), at which our results indicate small 
pressure effects on the fractionation (see Fig. |6ll, 
Wunder et al.l(l201 ih report a measurement of Li isotope 
fractionation between spodumene and aqueous fluid at 
r ~ 900K and P = 8 GPa to be -)-0.75 + 0.5 %c lower 
than the values measured at the same te mperature but 



lower pressures for the same systems in IWunder et al. 



(120061) . In order to check if we are able to reproduce 
this behavior with our method we computed the /? fac- 



tor of spodumene at high P = 8 GPa by using the lat- 
tice const ants of compressed spodumene determined by 



Arlt & Angel (2000 ). Because of the high compression, 
the resulting yS factor is 3 %c higher than the one de- 
rived for uncompressed solid. The /3 factor of fluid at 
the same {P, T) conditions increases by 1 .9 %c. This re- 
sults in pressure-driven decrease of the spodumene-fluid 
fractionation factor (A^Lispd._fluid) by 1.1 %o, w hich is in 



good agreement with the result of IWunder et al.l (1201 11) . 



4. 3. Fluid-mineral fractionation 

The different experiments on Li isotope fractiona- 
tion between Li -b earing minera l s and aqueou s solution 
at high P and T JWunder et all l2006i l2007h show the 
strongest enrichment in ^Li for staurolite and subse- 
quently lighter isotopic signatures for the fluid, mica 
and spodumene. An important test for our proposed 
computational method is to reproduce the sequence 
of experimentally observed fractionation factors. The 
crystal structures used in the calculation and the proce- 
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dure used to compute the y6 factors are described in pre- 
vious sections, and the relevant fi factors were already 
discussed. The computed fluid-mineral fractionation 
factors, A^Liniinerai-fluid, between staurolite, spodumene 
and mica, and aqueous solution are give n in figure [TOl 
togeth er with the experimental va lues of IWunder et alj 
(I2OO7I) for mica and staurolite and lWunder et alj (120061) 
for spodumene respectively. The errors of the computed 
fractionation factors are given in the figure caption and 
are derived assuming uncertainty in the computed vibra- 
tional frequencies coming from using BLYP functional, 
which is discussed in section [3j The computed curves 
correctly predict the fractionation sequence. The heavy 
Li isotope preferentially fractionates into staurolite with 
respect to aqueous solution, whereas spodumene is en- 
riched in ^Li. The computations also reproduce the ex- 
perimental results for both minerals on the quantitative 
level within 1-1.5 %o, taking into account the uncer- 
tainties in the calculated fractionation factors, and our 
prediction for spodumene is ideal. In case of mica the 
picture is more complicated as it has four polytypes and 
more than one Li substitution site. In figure[TO]we plot- 
ted the average mica-fluid fractionation factors com- 
puted for different polytypes. The resulted solid-fluid 
fractionation is higher than the experimental values by 
~ 1 - 2 %c, depending on the polytype. Nevertheless, 
our results confirm that among the considered miner- 
als, the fractionation between mica and the fluid is the 
smallest and that on average the mica containing mix- 
ture of different polytypes should be slightly enriched 
with light isotope comparing with fluid. We note that 
as the measurements for mica were performed at lower 
temperature of ~ 650 K, the error in the calculated frac- 
tionation factor between mica and fluid is significant 
and ~ 0.6 %o. The straightforward comparison of our 
results for mica crystalline solid with the experimental 
data is als o complicated as diffe rent reported measured 
samples of Wunder et al.l (120071) contained different rel- 



ative abundances of different polytypes. We also found 
that of all the crystalline solids considered here mica is 
the most sensitive to the change in the lattice parameters 
and computational setup. For instance, while the fi fac- 
tors for other minerals and the fluid are well converged 
(within 0. 1 %c) using the force constants obtained with 
the plane wave energy cutoff of 100 Ryd, the resulted 
fi factors for mica with this setup are overestimated by 
~ 1 .5 %c and the converged values were obtained by ap- 
plying much higher cutoff of 140 Ryd. 

We notice that for staurolite and mica the solid-fluid 
fractionation factors are overestimated by ~ 0.5 - 1 .5 %c. 
However, the lattice parameters used in the calculations 
of the crystalline solids are the one measured at am- 



bient conditions. At high temperatures solids undergo 
thermal expansion, which should result in the lower- 
ing of the fi factors. Observing the deviation of the 
computed solid-fluid fractionation factors for staurolite 
and mica we attempted to check for the effect of the 
thermal expansion of the lattice parameters of mod- 
eled solids on the derived fractionation factors. Ac- 



cording to the crystal structure data of ICameron et al 



(Il973h the lattice constants of spodumene expand by 
~ 0.5 % at r ~ 1000 K. Having such a pronounced ef- 
fect, we recalculated the fi factors of spodumene ?AT - 
573 K, 723 K and 10 33 K using temperatur e-dependent 
lattice parameters of ICameron et alJ (Il973l) . We found 
that the thermal expansion of spodumene results in 
~ 0.4 %c decrease in the fi factors for all the consid- 
ered te mperatures. Similar reduction is observed for 
micas. Russell & Gug genheimI (119991) showed that for 
the phlogopite IM mica the lattice parameters increase 
by ~ 0.5 % at r = 650 K. Assuming that Li-bearing 
micas undergo similar expansion we computed the fi 
factors with the lattice parameters rescaled by -1-0.5 %. 
The resulted fi factors are ~ 0.7 %c smaller, which indi- 
cates that inclusion of thermal expansion effect lowers 
the computed mica-fluid fractionation curves by 0.7 %o. 
We also computed the fi factor for staurolite assum- 
ing the ~ 0.5 % increase of its l attice parameters at 
1000 K dHoUand and PowellllioTTI) . The resulted^ fac- 
tor decreases by 0.6 %c. In addition we notice that the 
measurements for staur ohte were performed at higher 



pressure P = 3.5 GPa (IWunder et al.L l2007h . At such 



elevated pressure the fi factor for fluid increases by 
~ 0.5 %o (Fig. |6]l leading to the further decrease of the 
staurolite-fluid fractionation factor by 0.5 %c. The solid- 
fluid isotope fractionation factors resulted by applying 
the derived shifts in fi factors are given in the right panel 
of figure[TO] It is clearly visible that the corrections due 
to the thermal expansion of crystalline solids and the 
high pressure in case of staurolite make the prediction 
more consistent with the measurements. We note that 
on the right panel we plotted the solid-fluid fractiona- 
tion curves only at the temperature range corresponding 
to the experiment as being interested in direct compari- 
son of the computed values with the experiments we ap- 
plied the constant thermal expansion and pressure cor- 
rection derived only at these temperatures. The respec- 
tive corrections for staurolite and mica at other temper- 
atures may be different. 

Beside the thermal expansion effect and the uncer- 
tainties and systematic errors resulting from choice of 
the DFT functional there are other effects that could po- 
tentially increase the uncertainties in the calculated frac- 
tionation factors. These additional effects could arise 
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from the usage of the experimental equation of state for 
fluid and lattice parameters for crystalline solids, and 
uncertainties in the crystalline lattice site occupations 
as in the case of micas. On the other hand, it can not 
be guaranteed with full confidence that the ex perimen- 
tal measurements of IWunder et alj (120061 l2007i) . which 
indicate complete isotopic exchange, reflect an equi - 
librium fractionation (see for instance iLi et alJ (1201 ih ) 
[J. Nevertheless, the good agreement between theoreti- 
cal prediction and experimental data on the Li isotope 
fractionation between complex Li-bearing minerals and 
aqueous fluid shows that the outlined method for com- 
puting the isotope fractionation of fluids and crystals is 
a powerful tool, which can be successfully applied for 
prediction of isotopic signatures of complex Earth ma- 
terials under extreme conditions. 

5. Conclusions 

We propose a computationally eflicient approach for 
computation of the fi and isotope fractionation factors 
for complex minerals and fluids at high temperatures 
and pressures. We demonstrated that in order to derive 
the reliable fi factors for either minerals or fluids at high 
T and P it is sufficient to know the force constants acting 
on the substituted isotope. This reduces significantly the 
computational time and allows for computations of iso- 
tope fractionation in complex materials containing even 
hundreds of atoms. In case of fluids we show that the 
widely used technique of representing aqueous solution 
as an ion-hydration-shell cluster is not sufficient to re- 
produce the isotope fractionation in aqueous solutions 
at elevated temperatures and pressures, when the dy- 
namical character of the hydration shell and the com- 
pression of the fluid have to be accounted for. This can 
be achieved by ab initio molecular dynamics simulation 
technique, which allows for direct access to the dynam- 
ical distributions of water (fluid) molecules around the 
considered ion and proper consideration of compression 
effects. The relevant isotope fractionation factors can be 
computed on a set of uncorrected snapshot configura- 
tions extracted from the molecular dynamics trajectory. 

We show ffiat in the case of Li in aqueous solution it 
is sufficient to compute the y6 factors from the molecu- 
lar dynamics simulations performed with a simulation 
cefl containing a small number of atoms, which further 
reduces the computational time needed to perform the 



' Althou gh the argume nts supporting the equihbrium fractionation 
as given in IWunder et al J 12006, 2007) are convincing and the good 
agreement between our predictions and the measurements corroborate 
that scenario. 



task. A system containing a single Li, a charge com- 
pensating anion and 8 H2O molecules was sufficient to 
obtain the accurate fi factors within the uncertainties of 
the ab initio method used in the calculations. 

We verify our approach by computing the Li isotopes 
fractionation factors between Li-bearing minerals and 
aqueous solutions and their comparison with the exper- 
imental data. The computed Li fractionation factors be- 
tween staurolite, spodumene, mica and aqueous solu- 
tions reproduce the experimental results on quantitative 
and qualitative levels. We show that ab initio calcula- 
tions are able to predict the correct sequence of isotopes 
fractionation between considered materials as observed 
in the experiment. The computed fractionation factors 
are within 1 %o in agreement with the measured values. 
We also found that the thermal expansion of the solids 
aff'ects the isotope fractionation process and its inclu- 
sion improves the agreement with the experimental data. 

Our study shows that ab initio computer simulations 
represent a powerful tool for prediction and understand- 
ing of equilibrium stable isotope fractionation processes 
between various phases including aqueous solutions at 
high pressures and temperatures. We expect that with 
the increasing power of computers and performance of 
the computational software these methods will be exten- 
sively applied to complement analytical techniques and 
to interpret measured isotopic signatures. 
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